* Replication for "Particularly in the WEIRD World: Gender Gaps in Attitudes towards Homosexuality"
* Achim Hildebrandt & Sebastian Jäckle


* Version 16.12.2022

version 15

*ssc inst labutil
*ssc install g538schemes
clear all

*cd "YOUR-PATH\replication"
cd "E:\Eigene Dokumente\Freiburg\Forschung\Gender_Gap Homonegativity\replication"

*******************************************************************************
***																			***
***						A) Generate dataset for analysis					***
***																			***
*******************************************************************************

* Temporarily save Macro Data in dta-format
	* muslim population
	import excel "data\muslim_population.xlsx", sheet("muslims_for stata") firstrow clear
	rename Country country
	save "dump\muslim_population.dta", replace

	* Worldbank GNI p. Cap.
	import excel "data\Worldbank_GNI_p_Cap.xlsx", sheet("Worldbank_GNIpCap") firstrow clear
	rename CountryName country
	save "dump\GNI_p_Cap.dta", replace
	
	* UNDP HDI Human Development Index
	import excel "data\HDI_HDR2020_040722.xlsx", sheet("HDI_HDR2020_040722") firstrow clear
	replace country = "Bolivia" if country == "Bolivia (Plurinational State of)"
	replace country = "Hong Kong SAR" if country == "Hong Kong, China (SAR)"
	replace country = "South Korea" if country == "Korea (Republic of)"
	replace country = "Iran" if country == "Iran (Islamic Republic of)"
	replace country = "Russia" if country == "Russian Federation"
	replace country = "Vietnam" if country == "Viet Nam"
	replace country = "Venezuela" if country == "Venezuela (Bolivarian Republic of)"
	keep country hdi_2015-hdi_2019
	save "dump\HDI.dta", replace
	
	* UNDP GII Gender Inequality Index GII
	import excel "data\GII_HDR2020_040722.xlsx", sheet("GII_HDR2020_040722") firstrow clear
	replace country = "Bolivia" if country == "Bolivia (Plurinational State of)"
	replace country = "Hong Kong SAR" if country == "Hong Kong, China (SAR)"
	replace country = "South Korea" if country == "Korea (Republic of)"
	replace country = "Iran" if country == "Iran (Islamic Republic of)"
	replace country = "Russia" if country == "Russian Federation"
	replace country = "Vietnam" if country == "Viet Nam"
	replace country = "Venezuela" if country == "Venezuela (Bolivarian Republic of)"
	keep country gii_2015-gii_2019
	save "dump\GII.dta", replace
	
	* VDem Gender-Index
	import excel "data\Vdem_gender.xlsx", sheet("Vdem_gender") firstrow clear
	save "dump\Vdem_gender.dta", replace

	
	* Further Macro-Variables not tested in the article
	* WEF Gender Gap Index
	import excel "data\WEF Global_Gender_Gap_Index.xlsx", sheet("data") firstrow clear
	keep if Indicator == "Overall Global Gender Gap Index" & SubindicatorType == "Index"
	rename CountryName country
	save "dump\WEF_GGI.dta", replace

	* Social Watch Gender Equity Index (GEI) 2012: not continued after 2012
	import excel "data\Social Watch_GEI_2012.xlsx", sheet("GEI_2012_for_stata") firstrow clear
	rename Country country
	save "dump\SocialWatch_GEI.dta", replace
	
	* Worldbank Female Labor Force
	import excel "data\Worldbank_Female_Labor_Force.xlsx", sheet("Worldbank_Female_Labor_Force") firstrow clear
	rename Country country
	save "dump\FLF.dta", replace
	
	* UNDP GDI Gender Development Index
	import excel "data\GDI_HDR2020_040722.xlsx", sheet("GDI_HDR2020_040722") firstrow clear
	replace country = "Bolivia" if country == "Bolivia (Plurinational State of)"
	replace country = "Hong Kong SAR" if country == "Hong Kong, China (SAR)"
	replace country = "South Korea" if country == "Korea (Republic of)"
	replace country = "Iran" if country == "Iran (Islamic Republic of)"
	replace country = "Russia" if country == "Russian Federation"
	replace country = "Vietnam" if country == "Viet Nam"
	replace country = "Venezuela" if country == "Venezuela (Bolivarian Republic of)"
	keep country gdi_2015-gdi_2019
	save "dump\GDI.dta", replace
	
* Open WVS data
	use "data\WVS_Cross-National_Wave_7_stata_v4_0.dta"
	decode B_COUNTRY, gen(country)
	gen countrycode = B_COUNTRY_ALPHA

* Merge WVS with Macro Data
	merge m:1 country using "dump\muslim_population.dta", gen(_merge_muslim)
	merge m:1 country using "dump\GNI_p_Cap.dta", gen(_merge_GNI_p_Cap)
	merge m:1 country using "dump\HDI.dta", gen(_merge_HDI)
	merge m:1 country using "dump\GII.dta", gen(_merge_GII)
	merge m:1 country using "dump\Vdem_gender.dta", gen(_merge_VDem_gender)
	
	* Further Macro-Variables not tested in the article
	merge m:1 country using "dump\WEF_GGI.dta", gen(_merge_GGI)
	merge m:1 country using "dump\SocialWatch_GEI.dta", gen(_merge_GEI)
	merge m:1 country using "dump\FLF.dta", gen(_merge_FLF)
	merge m:1 country using "dump\GDI.dta", gen(_merge_GDI)
	drop if version == ""

* Missings for Macro Variables 
	* (see table_missing.xls for an overview)
	tab country if hdi_2018 ==.
	tab country if GNI_p_Cap_2018 ==.
	tab country if Vdem_Gender_2018 ==.
	tab country if PEW_muslim ==.
	tab country if gii_2018 ==.
	
	tab country if gdi_2018 ==.
	tab country if FLF_2018 ==.
	tab country if GGI_2018 ==.
	tab country if GEI ==.

* Recode and rename variables
	*dependent variable (homosexuality justifiable)
	rename Q182 homosexuality

	* religion
	recode Q289 (0=0) (5=1) (1 2 3 4 6 7 8 9 = 2), gen(religion)
	
	* Libya (no Q289 available --> recode Q289CS
	replace religion = 0 if country =="Libya" & Q289CS == 100000020
	replace religion = 1 if country =="Libya" & Q289CS < 100000000
	label define religion 0 "none" 1 "muslim" 2 "other religion"
	label values religion religion
	
	* importance of religion 
	recode Q6 (1=4) (2=3) (3=2) (4=1), gen(importance_religion)
	
	* HDI groups
	recode hdi_2018 (0.0001/0.549 = 1) (0.550/0.699 =2) (0.7/0.799 = 3) (0.8/1.0 = 4), gen(hdi_cat_2018)
	label define hdi_cat 1 "low" 2 "medium" 3 "high" 4 "very high"
	label values hdi_cat_2018 hdi_cat
	
	* Vdem-Gender groups
	recode Vdem_Gender_2018  (0.0000001/0.688 = 1) (0.689/0.798 = 2) (0.799/0.873 = 3) (0.874/1.0 = 4), gen(Vdem_Gender_cat_2018)
	label define Vdem_Gender_cat_2018 1 "Q1" 2 "Q2" 3 "Q3" 4 "Q4"
	label values Vdem_Gender_cat_2018 Vdem_Gender_cat_2018
	
	* male
	recode Q260 (1=1) (2=0), gen(male)
	label define male 1 "male" 0 "female"
	label values male male
	
	* age
	rename Q262 age
	
	* education - ISCED classification
	rename Q275 edu
	
	* Gender Role Attitudes (country mean)
	* Mean (Men make better political leaders than women do + University is more important for a boy than for a girl + Men make better business executives than women do)
	gen gender_attitude = (Q29 + Q30 + Q31)/3
	bysort country: egen mean_gender_attitude = mean(gender_attitude)
	
	* filter for all countries except for Andorra, Puerto Rico and Macau for which no Vdem-data is available 
	gen filter_country = 1
	replace filter_country = 0 if country == "Andorra"
	replace filter_country = 0 if country == "Puerto Rico"
	replace filter_country = 0 if country == "Macau SAR"
	
	* filter for all cases without missings for male age i.edu i.religion c.importance_religion
	gen filter_ind = 1 
	replace filter_ind = 0 if male == .
	replace filter_ind = 0 if age == .
	replace filter_ind = 0 if edu == .
	replace filter_ind = 0 if religion == .
	replace filter_ind = 0 if importance_religion == .

* label variables
	label variable homosexuality "Homosexuality justifiable"
	label variable male "Male"
	label variable age "Age"
	label variable edu "Highest educational level"
	label variable religion "Denomination"
	label variable importance_religion "Religiosity"
	label variable hdi_2018 "Human Development Index"
	label variable Vdem_Gender_2018 "Women Political Empowerment Index"
	label variable GNI_p_Cap_2018 "Gross National Income p. Cap (in US$)"
	label variable gii_2018 "Gender Inequality Index"
	label variable mean_gender_attitude "Mean gender role attitudes"
	label variable PEW_muslim "Muslim population (in %)"

save "dump\final_analysis_data.dta", replace

collapse (mean) gii = gii_2018 hdi = hdi_2018 GNI = GNI_p_Cap_2018 VDem = Vdem_Gender_2018 gender_att = mean_gender_attitude , by(countrycode)
save "dump\country_level_data.dta", replace

use "dump\final_analysis_data.dta"
collapse (mean) homo = homosexuality if filter_country==1, by(countrycode)
drop if countrycode == "EGY"
drop if countrycode == "TJK"
save "dump\mean_homo.dta", replace

 
	
*******************************************************************************
***																			***
***						B) Analysis											***
***																			***
*******************************************************************************

use "dump\final_analysis_data.dta"
set scheme plotplain
* Table A1: descriptive overview of all variables
estpost summarize homosexuality hdi_2018 Vdem_Gender_2018 male age  importance_religion  PEW_muslim  GNI_p_Cap_2018 mean_gender_attitude   if filter_country ==1 & filter_ind == 1, listwise  
esttab using "tables\table_A1a.rtf", cells("mean(fmt(2)) sd(fmt(2)) min(fmt(2)) max(fmt(2))") nomtitle  label replace

	* for those with missings
	estpost summarize homosexuality gii_2018 if filter_country ==1 & filter_ind == 1, listwise
	esttab using "tables\table_A1b.rtf", cells("mean(fmt(2)) sd(fmt(2)) min(fmt(2)) max(fmt(2))") nomtitle  label replace

	* for categorical variables
	tab religion if filter_country ==1 & filter_ind == 1
	tab edu if filter_country ==1 & filter_ind == 1
	
	* Figure 1
		* description of the dependent variable (homosexuality justifiable - means)

	* Means for Pakistan
	sum homosexuality if country == "Pakistan" & male == 1
	sum homosexuality if country == "Pakistan" & male == 0

	graph hbar homos [aweight = W_WEIGHT] if homos != . & filter_country == 1  , exclude0 ysc(r(1(1)10)) yla(1(1)10) over(country, sort(1) descending) xsize(4) ysize(7) ytitle("Homosexuality justifiable (mean)") blabel(total, format(%8.2f) size(tiny) position(base)) /// title("With weights")
	
	graph export "figures\Fig_01.emf", replace
	graph export "figures\Fig_01.pdf", replace
	
	* Figure 5
		* homosexuality by level of development and gender
			graph bar (mean) homosexuality if filter_country == 1  , exclude0  yscale(range(1 10))  over(male) over(hdi_cat_2018) asyvars blabel(total, format(%8.2f) size(vsmall)) legend(ring(0) position(10) bmargin(small) col(2) ) ytitle("Homosexuality justifiable (mean)") b1title("Human Development Index Group (UNDP)") bar(1, color(gs14)) bar(2, color(black))
			graph save "dump\Fig_05a.gph", replace
			
		* homosexuality by level of gender equality and gender
			graph bar (mean) homosexuality if filter_country == 1  , exclude0  yscale(range(1 10))  over(male) over(Vdem_Gender_cat_2018) asyvars blabel(total, format(%8.2f) size(vsmall))  legend(off) ytitle("") b1title("Gender Political Empowerment Index Quartile (V-Dem)") bar(1, color(gs14)) bar(2, color(black))
			graph save "dump\Fig_05b.gph", replace
			graph combine "dump\Fig_05a.gph" "dump\Fig_05b.gph", ycommo
			// scaleaxis edits
			gr_edit .plotregion1.graph1.scaleaxis.reset_rule 1 10 1 , tickset(major) ruletype(range) 
			gr_edit .plotregion1.graph2.scaleaxis.reset_rule 1 10 1 , tickset(major) ruletype(range) 
		
			graph export "figures\Fig_05.emf", replace
			graph export "figures\Fig_05.pdf", replace

	* Figure 3 		
		* Means of homosexuality justifiable for men and women			
	gen homosex_male = homosexuality if male == 1
	gen homosex_female = homosexuality if male == 0
	set scheme plotplain
	graph dot (mean) homosex_male (mean) homosex_female if filter_country ==1 & filter_ind == 1 & country != "Tajikistan" & country != "Egypt", exclude0 ysc(r(1(1)10)) yla(1(1)10)  scale(0.7) over(country, sort(2) descending) legend(order(1 "Male" 2 "Female") position(6) cols(2)) ysize(10) ytitle("Homosexuality justifiable (mean)") 
	graph export "figures\Fig_03.emf", replace		
	graph export "figures\Fig_03.pdf", replace
	
* Means of Homosexuality Justifiable for West & Non-Western countries
collapse homosex_male homosex_female if filter_country==1, by(countrycode)
gen west = 0
replace west = 1 if countrycode == "AUS"
replace west = 1 if countrycode == "CAN"
replace west = 1 if countrycode == "NZL"
replace west = 1 if countrycode == "USA"
replace west = 1 if countrycode == "DEU"
replace west = 1 if countrycode == "CYP"
replace west = 1 if countrycode == "GRC"
replace west = 1 if countrycode == "NLD"
replace west = 1 if countrycode == "ROU"
replace west = 1 if countrycode == "SRB"
sum homosex_male if west == 1
sum homosex_male if west == 0
sum homosex_female if west == 1
sum homosex_female if west == 0

by country: sum homosex_male

clear all
use "dump\final_analysis_data.dta"
	
* description of macro-level variables (not presented in the article)

	* HDI: UNDP 2018
	graph hbar hdi_2018 if filter_country == 1 , over(country, sort(1) descending) xsize(4) ysize(7) ytitle("HDI") blabel(total, format(%8.2f) size(tiny) position(base))
	graph export "additional_figures\Fig_1a.emf", replace
	graph export "additional_figures\Fig_1a.pdf", replace
	
	* Vdem_Gender: Women political empowerment index 2018
	graph hbar Vdem_Gender_2018 if filter_country == 1 , over(country, sort(1) descending) xsize(4) ysize(7) ytitle("Women Political Empowerment Index") 
	graph export "additional_figures\Fig_1b.emf", replace
	graph export "additional_figures\Fig_1b.pdf", replace

	* Muslim Population: PEW data for 2009
	graph hbar PEW_muslim if filter_country == 1 , over(country, sort(1) descending) xsize(4) ysize(7) ytitle("Muslim Population (in %)") ///title("PEW data")
	
	graph export "additional_figures\Fig_1c.emf", replace
	graph export "additional_figures\Fig_1c.pdf", replace
	
	* GNI per capita: Worldbank 2018
	graph hbar GNI_p_Cap_2018 if filter_country == 1 , over(country, sort(1) descending) xsize(4) ysize(7) ytitle("Gross National Income in current US$") 
	graph export "additional_figures\Fig_1d.emf", replace
	graph export "additional_figures\Fig_1d.pdf", replace

	* GII: UNDP 2018
	graph hbar gii_2018 if filter_country == 1 , over(country, sort(1) descending) xsize(4) ysize(7) ytitle("Gender Inequality Index") 
	graph export "additional_figures\Fig_1e.emf", replace
	graph export "additional_figures\Fig_1e.pdf", replace		

	*mean_gender_attitude
	graph hbar mean_gender_attitude if filter_country == 1 , over(country, sort(1) descending) xsize(4) ysize(7) ytitle("Index Mean Gender Role Attitudes") 
	graph export "additional_figures\Fig_1f.emf", replace
	graph export "additional_figures\Fig_1f.pdf", replace

		*further operationalizations not presented in the article
			* Female Labor Force: Worldbank 2018
			graph hbar FLF_2018 if filter_country == 1 , over(country, sort(1) descending) xsize(4) ysize(7) ytitle("Female Labor Force (Worldbank)") title("Female Labor Force 2018 (Worldbank)")
			graph export "additional_figures\Fig_1g.emf", replace
			graph export "additional_figures\Fig_1g.pdf", replace
			
			* Gender Gap Index: World Economic Forum 2018
			graph hbar GGI_2018 if filter_country == 1 , over(country, sort(1) descending) xsize(4) ysize(7) ytitle("Gender Gap Index (WEF)") title("Gender Gap Index 2018 (WEF)")
			graph export "additional_figures\Fig_1h.emf", replace
			graph export "additional_figures\Fig_1h.pdf", replace
			
			* GDI: UNDP 2018
			graph hbar gdi_2018 if filter_country == 1 , over(country, sort(1) descending) xsize(4) ysize(7) ytitle("Gender Development Index (UNDP)") title("Gender Development Index 2018 (UNDP)")
			graph export "additional_figures\Fig_1i.emf", replace
			graph export "additional_figures\Fig_1i.pdf", replace
			
			* Gender Equity Index: Social Watch 2012 (no later data available)
			graph hbar GEI if filter_country == 1 , over(country, sort(1) descending) xsize(4) ysize(7) ytitle("Gender Equity Index (Social Watch)") title("Gender Equity Index 2012 (Social Watch)")
			graph export "additional_figures\Fig_1j.emf", replace
			graph export "additional_figures\Fig_1j.pdf", replace
						
			* Muslim Population: Alternativ data from worldpopulationreview data for 2010
			graph hbar muslims_WPR, over(country, sort(1) descending) xsize(4) ysize(7) ytitle("Muslim Population (in %)")  title("World Population Review data")
			graph export "additional_figures\Fig_1k.emf", replace
			graph export "additional_figures\Fig_1k.pdf", replace

	
* Effect size
	
	*Cohens d
	sort country
	by country: esize  twosample homosexuality if filter_country == 1 , by(male) cohensd level(99)
	* save Cohen's d + ci lower bound & ci upper bound in new data set 
	statsby d=r(d) lb_d=r(lb_d) ub_d=r(ub_d) if filter_country == 1  & homosex != . , by(country) saving("dump\cohens_d.dta", replace): esize  twosample homosexuality, by(male) cohensd level(99) 
	
	*Cohens d (for correlation with homosexuality justifiable)
	sort countrycode
	by countrycode: esize  twosample homosexuality if filter_country == 1 , by(male) cohensd level(99)
	* save Cohen's d + ci lower bound & ci upper bound in new data set 
	statsby d=r(d) lb_d=r(lb_d) ub_d=r(ub_d) if filter_country == 1  & homosex != . , by(countrycode) saving("dump\cohens_d2.dta", replace): esize  twosample homosexuality, by(male) cohensd level(99)
	

	*Hedges g 
	sort country
	by country: esize  twosample homosexuality if filter_country == 1 , by(male) hedgesg level(99)
	* save hedges g + ci lower bound & ci upper bound in new data set 
	statsby g=r(g) lb_g=r(lb_g) ub_g=r(ub_g) if filter_country == 1  & homosex != . , by(country) saving("dump\hedges_g.dta", replace): esize  twosample homosexuality, by(male) hedgesg level(99)  
	
	* for construction of Hedges g/Cohens d graphs see lines 765-798


	
* Multilevel Models
	* M0: empty model
	mixed homosexuality || country: if filter_country ==1 & filter_ind == 1, variance reml level(99.9)
	eststo M0
	matrix A = e(N_g)
	estadd scalar ngrps = A[1,1]
	estat icc
	estadd scalar ICC = r(icc2)
	
	* M1: only individual level variables
	mixed homosexuality i.male age i.edu i.religion importance_religion || country: if filter_country ==1 & filter_ind == 1, variance reml level(99.9)
	eststo M1
	matrix A = e(N_g)
	estadd scalar ngrps = A[1,1]
	estat icc
	estadd scalar ICC = r(icc2)
	
	* M2: only individual level variables - male as random slope
	mixed homosexuality i.male age i.edu i.religion c.importance_religion || country: male if filter_country ==1 & filter_ind == 1 , variance reml level(99.9)
	eststo M2
	matrix A = e(N_g)
	estadd scalar ngrps = A[1,1]
	estat icc
	estadd scalar ICC = r(icc2)
	
	* M3: individual level + %muslim + hdi_2018 + Vdem_Gender_2018
	mixed homosexuality i.male age i.edu i.religion c.importance_religion PEW_muslim hdi_2018 Vdem_Gender_2018 || country: male if filter_country ==1 & filter_ind == 1, variance reml level(99.9)
	eststo M3
	matrix A = e(N_g)
	estadd scalar ngrps = A[1,1]
	estat icc
	estadd scalar ICC = r(icc2)	

	* M4 +Interaction male x %muslims &  male x hdi_2018  & male* Vdem_Gender_2018		
	use "dump\final_analysis_data.dta", replace
	mixed homosexuality i.male##c.PEW_muslim i.male##c.Vdem_Gender_2018 i.male##c.hdi_2018 age i.edu i.religion c.importance_religion  || country: male if filter_country ==1 & filter_ind == 1, variance reml level(99.9)
	eststo M4
	matrix A = e(N_g)
	estadd scalar ngrps = A[1,1]
	estat icc
	estadd scalar ICC = r(icc2)	
		
	* M5 +Interaction male x %muslims & triple Interaction male x hdi_2018 * Vdem_Gender_2018		
	use "dump\final_analysis_data.dta", replace
	mixed homosexuality i.male##c.PEW_muslim i.male##c.Vdem_Gender_2018##c.hdi_2018 age i.edu i.religion c.importance_religion  || country: male if filter_country ==1 & filter_ind == 1, variance reml level(99.9)
	eststo M5
	matrix A = e(N_g)
	estadd scalar ngrps = A[1,1]
	estat icc
	estadd scalar ICC = r(icc2)
	
*Table 4
	esttab M0 M1 M2 M3 M4 M5  using "tables\table_4.rtf", 	label nonumbers mgroups("M0" "M1" "M2" "M3" "M4" "M5", pattern( 1 1 1 1 1 1)) ///
															mtitles("empty model"   "only indiv.level" "+male as random slope" "+country level" "+interactions" "+ triple interaction male x HDI x VdemGender") ///
															interact(*) se aic scalars("ICC" "ngrps") ///
															transform(#*: exp(2*@) 2*exp(2*@)) equations(2:2:3:3:3:3, 3:3:4:4:4:4, .:.:2:2:2:2) stardrop(#*:) ///
															eqlabels(""  "Variance country-level" "Variance individual-level" "Variance Gender", none) replace 
															
															** --> the table additionally includes R2-Mikro and R2 Makro (Snijders/Bosker 1994)

	
	set scheme plotplain
* AME single graphs

	*Figure 6
	* Average Marginal Effect of male by HDI_2018
		margins, dydx(male) at (hdi_2018=(0.3(0.05)1)) 
		marginsplot, recast(line) recastci(rarea) ciopt(color(%20)) title("") xtitle("Human Development Index (UNDP)") ytitle("Average Marginal Effect of Gender (Male)")
		graph save "dump\Fig_06a.gph", replace
		
	* Average Marginal Effect of male by Vdem_Gender_2018
		margins, dydx(male) at (Vdem_Gender_2018=(0.3(0.05)1)) 
		marginsplot, recast(line) recastci(rarea) ciopt(color(%20)) title("") xtitle("Women Political Empowerment Index (V-Dem)") ytitle("")
		graph save "dump\Fig_06b.gph", replace
		
	* Average Marginal Effect of male by Muslim Population
		margins, dydx(male) at (PEW_muslim=(0(10)100)) 
		marginsplot, recast(line) recastci(rarea) ciopt(color(%20)) title("") xtitle("Muslim population in %") ytitle("")
		graph save "dump\Fig_06c.gph", replace

graph combine  "dump\Fig_06a.gph" "dump\Fig_06b.gph" "dump\Fig_06c.gph", col(3) xsize(7) ysize(4) scale(0.8) ycommon
graph export "figures\Fig_06.emf", replace
graph export "figures\Fig_06.pdf", replace		

	* for robustness graphs only:
		* Average Marginal Effect of male by HDI_2018 (for robustness)
			margins, dydx(male) at (hdi_2018=(0.3(0.05)1)) 
			marginsplot, recast(line) recastci(rarea) ciopt(color(%20)) title("Main model M5" "with Human Development Index (UNDP") xtitle("Human Development Index (UNDP)") ytitle("Average Marginal Effect of Gender (Male)")
			graph save "dump\Fig_A1a.gph", replace
		
		* Average Marginal Effect of male by Vdem_Gender_2018 (for robustness)
			margins, dydx(male) at (Vdem_Gender_2018=(0.3(0.05)1)) 
			marginsplot, recast(line) recastci(rarea) ciopt(color(%20)) title("Main model M5" "with Women Political Empowerment Index (V-Dem)") xtitle("Women Political Empowerment Index") ytitle("Average Marginal Effect of Gender (Male)")
			graph save "dump\Fig_A2a.gph", replace

	* Figure 7
	* AME Contourplot
		set matsize 10000
	
		margins, dydx(male) at(hdi_2018=(0.3(0.1)1) Vdem_Gender_2018=(0.3(0.1)1) ) vsquish saving("dump/predictions_M4.dta", replace)
		use "dump/predictions_M4.dta", clear
		merge 1:1 _n using "dump\country_level_data.dta"
		
		rename _at3 Vdem_Gender_2018
		rename _at4 hdi_2018
		rename _margin homosexuality
		
		set scheme plottigblind
		
		gen pos = 3
		replace pos = 9 if countrycode == "PER"
		replace pos = 9 if countrycode == "MNG"
		replace pos = 9 if countrycode == "TWN"
		replace pos = 9 if countrycode == "GRC"
		replace pos = 12 if countrycode == "CAN"
		replace pos = 9 if countrycode == "KOR"
		replace pos = 6 if countrycode == "USA"
		replace pos = 5 if countrycode == "HKG"
		replace pos = 1 if countrycode == "NLD"
		replace pos = 9 if countrycode == "MAR"
		replace pos = 9 if countrycode == "KGZ"
			
		twoway (contour homosexuality Vdem_Gender_2018 hdi_2018,  ccuts(-.8(.1).5) ytitle("Women Political Empowerment Index (V-Dem)") xtitle("Human Development Index (UNDP)") ztitle("Average Marginal Effect Gender (Male)")  xlab(0.3(.1)1) ylab(0.3(.1)1) )  || ///
		(scatter VDem hdi, mcolor(black) msymbol(dh) mlabel(countrycode) mlabsize(tiny) mlabcolor(black) mlabpos(3) mlabv(pos) mlabgap(*.5) ) ,xsize(10) ysize(6) scale(1) plotr(m(zero)) 
		graph export "figures\Fig_07_color.emf", replace
		graph export "figures\Fig_07_color.pdf", replace

		twoway (contour homosexuality Vdem_Gender_2018 hdi_2018, crule(linear) scolor(white) ecolor(black) ccuts(-.8(.1).5) ytitle("Women Political Empowerment Index (V-Dem)") xtitle("Human Development Index (UNDP)") ztitle("Average Marginal Effect Gender (Male)")  xlab(0.3(.1)1) ylab(0.3(.1)1) )  || ///
		(scatter VDem hdi, mcolor(black) msymbol(dh) mlabel(countrycode) mlabsize(tiny) mlabcolor(black) mlabpos(3) mlabv(pos) mlabgap(*.5) ) ,xsize(10) ysize(6) scale(1) plotr(m(zero)) 
		graph export "figures\Fig_07_grey.emf", replace
		graph export "figures\Fig_07_grey.pdf", replace
		
** Robustness checks with different operationalizations for HDI and Women Empowerment		

* Models for Table A5
	use "dump\final_analysis_data.dta", replace
	mixed homosexuality i.male##c.PEW_muslim male##c.Vdem_Gender_2018##c.GNI_p_Cap_2018 age i.edu i.religion c.importance_religion  || country: male if filter_country ==1 & filter_ind == 1, variance reml level(99.9)
	eststo M5a
	matrix A = e(N_g)
	estadd scalar ngrps = A[1,1]
	estat icc
	estadd scalar ICC = r(icc2)
	
	mixed homosexuality i.male##c.PEW_muslim male##c.mean_gender_attitude##c.hdi_2018 age i.edu i.religion c.importance_religion  || country: male if filter_country ==1 & filter_ind == 1, variance reml level(99.9)
	eststo M5b
	matrix A = e(N_g)
	estadd scalar ngrps = A[1,1]
	estat icc
	estadd scalar ICC = r(icc2)
	
	mixed homosexuality i.male##c.PEW_muslim male##c.gii_2018##c.GNI_p_Cap_2018 age i.edu i.religion c.importance_religion  || country: male if filter_country ==1 & filter_ind == 1, variance reml level(99.9)
	eststo M5c	
	matrix A = e(N_g)
	estadd scalar ngrps = A[1,1]
	estat icc
	estadd scalar ICC = r(icc2)

	esttab M5a M5b M5c using "tables\table_A1.rtf", 	label nonumbers mgroups("M5a" "M5b" "M5c", pattern( 1 1 1 )) ///
															mtitles("GNI p.Cap."   "Mean Gender Role Attitude Index" "Gender Inequality Index" ) ///
															interact(*) se aic scalars("ICC" "ngrps") ///
															transform(#*: exp(2*@) 2*exp(2*@))  equations(4:4:4, 2:2:2, 3:3:3 ) stardrop(#*:) ///
															eqlabels(""  "Variance Gender" "Variance country level" "Variance individual level"  , none) replace 
															
															** --> the table additionally includes R2-Mikro and R2 Makro (Snijders/Bosker 1994) see Excel-File "R2 mikro-makro.xlsx" for the calculation  

* Models for Figures A1-A5
* M5a +Interaction male x %muslims & triple Interaction male x GNI_p_Cap_2018 * Vdem_Gender_2018		
	use "dump\final_analysis_data.dta", replace
	mixed homosexuality i.male##c.PEW_muslim male##c.Vdem_Gender_2018##c.GNI_p_Cap_2018 age i.edu i.religion c.importance_religion  || country: male if filter_country ==1 & filter_ind == 1, variance reml level(99.9)
	eststo M5a

	set scheme plotplain
	* AME einzeln
	* Average Marginal Effect of male by GNI_p_Cap_2018
		margins, dydx(male) at (GNI_p_Cap_2018=(0(10000)80000)) 
		marginsplot, recast(line) recastci(rarea) ciopt(color(%20)) title("") xtitle("GNI p. Cap. (in US$)") ytitle("") title("Alternative Operationalization:" "GNI p. Cap. (Worldbank)")
		graph save "dump\Fig_A1b.gph", replace

	* AME Contourplot
		set matsize 10000
	
		margins, dydx(male) at(GNI_p_Cap_2018=(0(10000)80000) Vdem_Gender_2018=(0.3(0.1)1)) vsquish saving("dump/predictions_M5a.dta", replace)
		use "dump/predictions_M5a.dta", clear
		merge 1:1 _n using "dump\country_level_data.dta"
		
		rename _at3 Vdem_Gender_2018
		rename _at4 GNI_p_Cap_2018
		rename _margin homosexuality
		
		set scheme plottigblind
		
		gen pos = 3
		replace pos = 9 if countrycode == "PER"
		replace pos = 9 if countrycode == "MNG"
		replace pos = 9 if countrycode == "TWN"
		replace pos = 9 if countrycode == "GRC"
		replace pos = 12 if countrycode == "CAN"
		replace pos = 9 if countrycode == "KOR"
		replace pos = 6 if countrycode == "USA"
		replace pos = 5 if countrycode == "HKG"
		replace pos = 1 if countrycode == "NLD"
		replace pos = 5 if countrycode == "KGZ"
		replace pos = 10 if countrycode == "TUN"
		replace pos = 6 if countrycode == "PER"
		replace pos = 2 if countrycode == "VNM"
		replace pos = 11 if countrycode == "ZWE"
		replace pos = 9 if countrycode == "KEN"
		replace pos = 12 if countrycode == "ETH"
		replace pos = 4 if countrycode == "GTM"
		replace pos = 6 if countrycode == "NGA"
			
		twoway (contour homosexuality Vdem_Gender_2018 GNI_p_Cap_2018,  ccuts(-1.2(.2)2.2) ytitle("Women Political Empowerment Index (V-Dem)") xtitle("GNI p. Cap. (in US$)") ztitle("Average Marginal Effect of Gender (Male)")  xlab(0(10000)80000) ylab(0.3(.1)1) )  || ///
		(scatter VDem GNI, mcolor(black) msymbol(dh) mlabel(countrycode) mlabsize(tiny) mlabcolor(black) mlabpos(3) mlabv(pos) mlabgap(*.5) ) ,xsize(10) ysize(6) scale(1) plotr(m(zero))
		graph export "figures\Fig_A3.emf", replace
		graph export "figures\Fig_A3.pdf", replace

	
	
* M5b +Interaction male x %muslims & triple Interaction male x hdi_2018 * mean_gender_attitude		
	use "dump\final_analysis_data.dta", replace
	mixed homosexuality i.male##c.PEW_muslim male##c.mean_gender_attitude##c.hdi_2018 age i.edu i.religion c.importance_religion  || country: male if filter_country ==1 & filter_ind == 1, variance reml level(99.9)
	eststo M5b
	set scheme plotplain
	* AME einzeln
		
	* Average Marginal Effect of male by mean_gender_attitude
		margins, dydx(male) at (mean_gender_attitude=(1(0.2)4)) 
		marginsplot, recast(line) recastci(rarea) ciopt(color(%20)) title("") xtitle("Mean Gender Role Attitudes") ytitle("") title("Alternative Operationalization:" "Mean Gender Role Attitudes Index (WVS)")
		graph save "dump\Fig_A2b.gph", replace

	* AME Contourplot
		set matsize 10000
	
		margins, dydx(male) at(hdi_2018=(0.3(0.1)1) mean_gender_attitude=(1(0.2)4)) vsquish saving("dump/predictions_M5b.dta", replace)
		use "dump/predictions_M5b.dta", clear
		merge 1:1 _n using "dump\country_level_data.dta"
		
		rename _at3 mean_gender_attitude
		rename _at4 hdi_2018
		rename _margin homosexuality
		
		set scheme plottigblind
		
		gen pos = 3
		replace pos = 12 if countrycode == "PER"
		replace pos = 9 if countrycode == "MNG"
		replace pos = 9 if countrycode == "TWN"
		replace pos = 9 if countrycode == "GRC"
		replace pos = 3 if countrycode == "CAN"
		replace pos = 9 if countrycode == "KOR"
		replace pos = 6 if countrycode == "USA"
		replace pos = 5 if countrycode == "HKG"
		replace pos = 1 if countrycode == "NLD"
		replace pos = 9 if countrycode == "MAR"
		replace pos = 3 if countrycode == "KGZ"
		replace pos = 11 if countrycode == "NZL"
		replace pos = 6 if countrycode == "CHN"
		replace pos = 12 if countrycode == "IDN"
		replace pos = 9 if countrycode == "MYS"
		replace pos = 6 if countrycode == "LBN"
		replace pos = 9 if countrycode == "BRA"
			
		twoway (contour homosexuality mean_gender_attitude hdi_2018,  ccuts(-1.1(.1).5) ytitle("Mean Gender Role Attitudes (WVS)") xtitle("Human Development Index (UNDP)") ztitle("Average Marginal Effect of Gender (Male)")  xlab(0.3(.1)1) ylab(1(.25)4) )  || ///
		(scatter gender_att hdi, mcolor(black) msymbol(dh) mlabel(countrycode) mlabsize(tiny) mlabcolor(black) mlabpos(3) mlabv(pos) mlabgap(*.5) ) ,xsize(10) ysize(6) scale(1) plotr(m(zero))
		graph export "figures\Fig_A4.emf", replace
		graph export "figures\Fig_A4.pdf", replace	

* M5c +Interaction male x %muslims & triple Interaction male x GNI_p_Cap_2018 * gii_2018		
	use "dump\final_analysis_data.dta", replace
	mixed homosexuality i.male##c.PEW_muslim male##c.gii_2018##c.GNI_p_Cap_2018 age i.edu i.religion c.importance_religion  || country: male if filter_country ==1 & filter_ind == 1, variance reml level(99.9)
	eststo M5c
	set scheme plotplain	
	
	* Average Marginal Effect of male by gii_2018
		margins, dydx(male) at (gii_2018=(0(0.02)0.6)) 
		marginsplot, recast(line) recastci(rarea) ciopt(color(%20)) title("Alternative Operationalization:" "Gender Inequality Index (UNDP)") xtitle("Gender Inequality Index") ytitle("")
		graph save "dump\Fig_A2c.gph", replace
	
	* AME Contourplot
		set matsize 10000
	
		margins, dydx(male) at(GNI_p_Cap_2018=(0(10000)80000) gii_2018=(0(0.02)0.6) ) vsquish saving("dump\predictions_M5c.dta", replace)
		use "dump\predictions_M5c.dta", clear
		merge 1:1 _n using "dump\country_level_data.dta"
		
		rename _at3 gii_2018
		rename _at4 GNI_p_Cap_2018
		rename _margin homosexuality
		
		set scheme plottigblind
		
		gen pos = 3
		replace pos = 9 if countrycode == "PER"
		replace pos = 9 if countrycode == "NIC"
		replace pos = 9 if countrycode == "EGY"
		replace pos = 4 if countrycode == "JOR"
		replace pos = 9 if countrycode == "MMR"
		replace pos = 6 if countrycode == "ETH"
		replace pos = 9 if countrycode == "PAK"
		replace pos = 9 if countrycode == "LBN"
		replace pos = 9 if countrycode == "TWN"
		replace pos = 9 if countrycode == "GRC"
		replace pos = 3 if countrycode == "CAN"
		replace pos = 9 if countrycode == "KOR"
		replace pos = 5 if countrycode == "HKG"
		replace pos = 10 if countrycode == "MAR"
		replace pos = 3 if countrycode == "KGZ"
		replace pos = 12 if countrycode == "IDN"
		replace pos = 9 if countrycode == "MYS"
		replace pos = 6 if countrycode == "CAN"
		replace pos = 9 if countrycode == "UKR"
		replace pos = 7 if countrycode == "BOL"
			
		twoway (contour homosexuality gii_2018 GNI_p_Cap_2018,  ccuts(-2.3(.2).1) ytitle("Gender Inequality Index (UNDP)") xtitle("GNI p. Cap. (in US$)") ztitle("Average Marginal Effect of Gender (Male)")  xlab(0(10000)80000) ylab(0(.05)0.6) )  || ///
		(scatter gii GNI, mcolor(black) msymbol(dh) mlabel(countrycode) mlabsize(tiny) mlabcolor(black) mlabpos(3) mlabv(pos) mlabgap(*.5) ) ,xsize(10) ysize(6) scale(1) plotr(m(zero))
		graph export "figures\Fig_A5.emf", replace
		graph export "figures\Fig_A5.pdf", replace	

* graphs for robustness checks /alternative operationalizations of development & gender equality

	* Fig A1
	graph combine  "dump\Fig_A1a.gph" "dump\Fig_A1b.gph" , col(2) xsize(7) ysize(4) scale(0.8) ycommon
	graph export "figures\Fig_A1.emf", replace
	graph export "figures\Fig_A1.pdf", replace

	* Fig A2
	graph combine  "dump\Fig_A2a.gph" "dump\Fig_A2b.gph" "dump\Fig_A2c.gph"  , col(3) xsize(7) ysize(4) scale(0.8) ycommon
	graph export "figures\Fig_A2.emf", replace
	graph export "figures\Fig_A2.pdf", replace




	
* Effect size graph

* Cohens d
clear all
use "dump\cohens_d.dta"		
set scheme plotplain

	* Figure 2
	* Dotplot + CI (reversed Cohens D for male - female)
	gen d2 = d*-1
	gen lb2_d =lb_d*-1 
	gen ub2_d =ub_d*-1
	sort d2
	gen country3 = _n
	sort country
    labmask country3, val(country)
	twoway (rspike  ub2_d lb2_d country3, legend(off) horizontal xline(0, lcolor(red))) (scatter country3 d2) , scale(0.7) ylabel(1(1)54, valuelabel) legend( order(1)) ytitle("") xsize(4) ysize(6) xtitle("Cohen's d {&plusmn} 99% CI")
	graph export "figures\Fig_02.emf", replace
	graph export "figures\Fig_02.pdf", replace
	
* Correlation Cohen's d & mean homosexuality justifiable
	clear all
	use "dump\cohens_d2.dta"	
	merge 1:1 countrycode using "dump\mean_homo.dta"
	gen d2 = d*-1
	corr d2 homo
	
	* Figure 4
	* Cohen's d vs. mean homosexuality justifiable

		gen pos = 3
		replace pos = 10 if countrycode == "SRB"
		replace pos = 1 if countrycode == "TWN"
		replace pos = 6 if countrycode == "ECU"
		replace pos = 11 if countrycode == "THA"
		replace pos = 7 if countrycode == "VEN"
		replace pos = 1 if countrycode == "MYS"
		replace pos = 4 if countrycode == "KOR"
		replace pos = 12 if countrycode == "KEN"
		replace pos = 9 if countrycode == "ARM"
		replace pos = 6 if countrycode == "TUN"
		replace pos = 9 if countrycode == "IDN"
		replace pos = 6 if countrycode == "KGZ"
		replace pos = 6 if countrycode == "LBY"
		replace pos = 1 if countrycode == "JOR"

scatter d2 homo, mcolor(black) msymbol(oh) mlabel(countrycode) mlabsize(tiny) mlabcolor(black) mlabpos(3) mlabv(pos) mlabgap(*.5)  ytitle("Effect size for gender: Cohen's d" "(male minus female)") xtitle("Mean of homosexuality justifiable") xlab(1(1)10)
gr_edit .xaxis1.edit_tick 10 1 `""1" " never" "justifiable""', tickset(major)
gr_edit .xaxis1.edit_tick 9 10 `""10" "always" "justifiable""', tickset(major)
graph export "figures\Fig_04.emf", replace
graph export "figures\Fig_04.pdf", replace

	


* Alternatively: Hedges g
clear all
use "dump\hedges_g.dta"		
set scheme plotplain
	* bar graph with significannt Hedges g in green color
	graph hbar (asis) g, over(country, sort(1)) bar(1, color(grey%20)) xsize(4) ysize(6) scale(0.7)	ytitle("Hedges g")	

	forvalues i = 1(1)17 {
	gr_edit .plotregion1.bars[`i'].EditCustomStyle , j(-1) style(shadestyle(color(dkgreen)))
	gr_edit .plotregion1.bars[`i'].EditCustomStyle , j(-1) style(shadestyle(color(%100)))
	gr_edit .plotregion1.bars[`i'].EditCustomStyle , j(-1) style(linestyle(color(dkgreen)))
	gr_edit .plotregion1.bars[`i'].EditCustomStyle , j(-1) style(linestyle(color(%100)))
	// bars[57] color
	}

	gr_edit .plotregion1.bars[54].EditCustomStyle , j(-1) style(shadestyle(color(dkgreen)))
	gr_edit .plotregion1.bars[54].EditCustomStyle , j(-1) style(shadestyle(color(%100)))
	gr_edit .plotregion1.bars[54].EditCustomStyle , j(-1) style(linestyle(color(dkgreen)))
	gr_edit .plotregion1.bars[54].EditCustomStyle , j(-1) style(linestyle(color(%100)))
	// bars[1] color	
	graph export "additional_figures\Fig_02_HedgesG_BAR.emf", replace	
	graph export "additional_figures\Fig_02_HedgesG_BAR.pdf", replace		

	
	* Dotplot + CI
	sort g
	gen country2 = _n
	sort country
    labmask country2, val(country)
	twoway (rspike lb_g ub_g country2, legend(off) horizontal xline(0, lcolor(red))) (scatter country2 g) , scale(0.7) ylabel(1(1)54, valuelabel) legend( order(1)) ytitle("") xsize(4) ysize(6) xtitle("Hedges g {&plusmn} 99% CI") 
	graph export "additional_figures\Fig_04_HedgesG.emf", replace
	graph export "additional_figures\Fig_04_HedgesG.pdf", replace
